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Abstract — To  a  great  extent,  the  success  of  advanced  image-guided  med¬ 
ical  procedures  hinges  on  non-rigid  volume  registration.  For  example,  non- 
rigid  registration  must  be  applied  in  interventional  approaches  where  intra¬ 
operative  information  is  used  to  update  high-quality  preoperative  data;  in 
follow-up  studies  in  order  to  assess  time-evolution  of  development;  aging, 
pathology  or  treatment;  and  in  many  other  applications  including  inter¬ 
subject  variability  and  population-based  atlas  construction. 

In  this  paper  we  examine  several  computational  schemes  that  warps  one 
volumetric  dataset  onto  another.  We  also  explore  the  inevitable  trade-off 
between  the  computational  load  and  the  incorporation  of  sophisticated  sim¬ 
ilarity  measures  necessary  for  multimodal  volumes.  Estimated  deforma¬ 
tion  fields  are  based  on  the  variational  formulation  of  Partial  Derivative 
Equations  (PDEs),  which  includes  a  similarity  and  a  regularization  term. 
We  compare  numerical  solutions  to  this  problem  using  the  Euler-Lagrange 
equations  (EL),  the  Finite  Elements  discretization  (FE),  and  Decoupled  Op¬ 
timization  over  the  possible  deformations  (DO). 

Keywords — Registration,  Correlation  Coefficient,  Finite  Elements  Method, 
Euler-Lagrange,  Template  Matching,  Optimization 

I.  INTRODUCTION 

Advanced  computer-assisted  interventions  require  3D  warp¬ 
ing  methods  to  perform  non-rigid  registrations  among  different 
volume  datasets,  an  essential  feature  of  interventional  proce¬ 
dures  [1]  where  preoperative  data  and  models  must  be  updated 
using  intraoperative  imaging  to  successfully  manage  natural  and 
induced  alterations  of  anatomy.  In  cases  such  as  these,  time  con¬ 
suming  image  acquisition  and  processing  is  curtailed  by  oper¬ 
ating  room  procedures,  and  this  obviously  limits  the  guidance 
that  intraoperative  images  can  provide.  In  an  effort  to  over¬ 
came  this  obstacle,  we  suggest  warping  high  quality  preoper¬ 
ative  data  and  models  into  the  coordinate  system  of  intraoper¬ 
ative  real  anatomy  (scanned  by  an  interventional  imaging  de¬ 
vice),  which,  in  turn,  is  used  to  guide  the  surgical  procedure. 
The  warp  is  obtained  by  estimating  a  deformation  field  and  then 
non-rigidly  adapting  the  proeoperative  data  to  the  interoperative 
data  by  non-rigidly  registering  the  preoperative  and  intraopera¬ 
tive  datasets.  Additional  applications  include  follow-up  studies 
of  the  same  individual  by  non-rigid  registration  of  image  data 
(e.g.,  developmental  [2]  and  aging  effects  [3]  studies,  degener¬ 
ative  pathology  [4]  and  treatment  [5]  evolution).  Nonrigid  reg¬ 
istration  is  also  crucial  to  intrasubject  comparison  (i.e.,  in  the 
construction  of  population-based  atlas  [6]  where  datasets  from 
different  individuals  must  be  related  to  a  canonical  frame).  Fur¬ 
ther,  non-rigid  registration  must  incorporate  a  priory  anatomic 
knowledge  into  automatic  segmentation  algorithms,  such  as  in 
the  framework  of  template-driven  segmentation  proposed  by 
Warfield  et  al.  [7]. 

Several  approaches  have  been  proposer  for  volume  data  regis¬ 
tration  that  rely  on  a  multi-resolution  decomposition  of  the  vol¬ 
umes.  With  regard  to  non-rigid  registration  algorithms,  the  for¬ 
mulas  must  recognize  potential  matches  based  on  the  similarity 
of  the  neighborhoods  to  corresponding  points  in  both  datasets. 
From  a  mathematical  standpoint,  these  algorithms  are  typically 
encoded  in  such  a  way  that  similarity  functions  between  neigh¬ 
borhoods  are  defined.  In  general,  many  combinations  are  pos¬ 
sible  for  acceptable  levels  of  the  similarity  function  ,  although 
the  inverse  problem  of  estimating  the  deformation  field  between 
both  datasets  is  ill-conceive,  a  typical  problem  in  regularization 
theory  as  defined.  For  example,  to  constrain  the  number  of  so¬ 


lutions  into  a  single  one  solution. 

An  optimal  way  of  including  both  the  similarity  and  the  reg¬ 
ularization  conditions  in  an  operating  equation  is  to  use  mod¬ 
els  based  on  variational  formulations  and  on  PDEs.  For  exam¬ 
ple,  the  well-known  sum  of  squared  differences  similarity  func¬ 
tion  is  used  in  [8]  with  a  viscous  fluid  constraint  that  does  not 
need  any  small  deformation  assumption,  in  [9]  with  an  elasticity 
constraint,  solving  the  PDE  with  the  FE  method.  Dengler  and 
Schmidt  [10]  proposed  an  optical  flow  method  with  an  elasticity 
constraint,  solved  with  a  multiresolution  pyramid  with  FEM  dis¬ 
cretization  (utilizing  the  sign  of  Laplacian  pyramid  as  the  match¬ 
ing  features).  Another  useful  similarity  measure  is  the  correla¬ 
tion  index,  which  considers  neighborhoods  to  be  similar  if  their 
intensities  are  related  by  an  affine  function.  This  measure  was 
applied  in  [11]  with  an  elasticity  constraint  as  regularization; 
solving  the  PDE  with  the  Jacobi  method;  unlikewise  in  [12]  but 
solving  the  PDE  iteratively  with  FE. 

Though  there  has  been  significant  efforts  to  solve  the  nonrigid 
registration  problem,  there  has  been  relatively  little  research  into 
more  complex  similarity  functions,  such  as  the  mutual  informa¬ 
tion  broadly  used  for  rigid  registration  [13],  which  are  essential 
for  multimodal  registration.  However,  arbitrary  similarity  func¬ 
tions  (see  [14]  for  examples)  are  difficult  to  be  incorporated  into 
the  computational  algorithms  due  to  their  lack  of  mathematical 
tractability.  This  deficiency,  fortunately,  can  be  compensated 
for  by  directly  computing  the  similarity  (as  opposed  to  apply¬ 
ing  a  an  analytical  derivation  of  the  similarity  function).  In  this 
paper,  we  will  frame  the  variational  problem  as  one  of  “equiva¬ 
lent  optimization”,  and  solve  it  by  determining  the  set  of  voxel 
displacements.  We  will  the  compare  these  results  with  those  ob¬ 
tained  using  the  Euler  Lagrange  equations  and  Finite  Elements 
discretizations. 

II.  VARIATIONAL  FORMULATION 
A  common  approach  to  nonrigid  registration  lies  in  establish¬ 
ing  the  displacement  field  d{f)  which  minimizes  an  energy  func¬ 
tional  with  a  similarity  and  a  regularization  term  relating  to  ref¬ 
erence  and  current  volumes;  that  is, 

J{u{f))  =  [  [S(f,u(r))  +  1Z{r,u{r))\dr  (1) 

where  S  is  the  similarity  measure  and  1Z  is  the  regularization 
term.  1Z  is  usually  a  smoothness  constraint. 

In  order  to  make  comparisons,  we  have  chosen  the  correla¬ 
tion  coefficient  as  the  similarity  function  and  the  sum  of  squared 
partial  derivatives  of  the  displacement  as  the  regularization  con¬ 
straint.  Hence,  the  correlation  coefficient  between  two  random 
variables  x  and  y  is  defined  as: 

CC(x,y)2  =  [^]2  (2) 

G y 

Where  axy ,  is  the  cross  covariance  and  ax ,  ay  are  the  standard 
deviations  of  x  and  y  respectively.  Inherit  in  the  non-rigid  reg¬ 
istration  problem,  are  the  intensities  of  two  tentatively  corre¬ 
sponding  voxels  x  and  y.  The  correlation  coefficient  is  esti¬ 
mated  from  two  windows  centered  on  those  voxels.  However, 
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A.  Euler-Lagrange 


window  size  must  be  traded  off  in  order  to  be  local  and  to  har¬ 
vest  a  sufficient  number  of  samples  for  a  meaningful  estimation. 
The  correlation  coefficient  facilitates  similarity  among  voxels  if 
the  intensities  inside  are  related  through  an  affine  function. 

For  the  regularization  term,  we  use 

R  =  |V  ®  u(f)ffro'  (3) 

where  u  is  the  deformation  field,  and  which  is  extensively 
used  in  optical  flow  applications  [15]  and  essentially  imposes 
a  smoothness  constraint  on  the  deformation  field. 

III.  Decoupled  Optimization 

From  a  numerical  point  of  view,  the  minimization  of  the  func¬ 
tional  (1),  i.e.,  the  estimation  of  the  displacement  field,  can  be 
viewed  as  an  optimization  problem,  since  the  displacement  vec¬ 
tors  for  every  voxel  are  the  parameters  that  must  be  optimized. 
Considered  in  this  way,  the  functional  becomes  a  function  or 
cost  of  the  set  of  displacements;  hence,  the  minimization  of  the 
functional  is  equivalent  to  the  minimization  of  a  multi-variate 
cost.  This  relationship,  in  turn,  allows  us  to  conduct  a  numeri¬ 
cal  search  for  the  minimum,  thus  avoiding  the  mathematical  in¬ 
tractability  of  complicated  functionals.  Nevertheless,  optimiza¬ 
tion  is  a  daunting  task  that  requires  a  clever  approach.  Notice  the 
number  of  unknown  parameters  is  three  times  (in  3D)  the  num¬ 
ber  of  voxels  and  the  highly  nonlinear  nature  of  the  functionals 
with  nontrivial  similarity  or  regularization  terms. 

In  this  paper,  we  propose  here  a  numerical  method,  general 
enough  to  deal  with  any  similarity  function  or  regularization 
term  based  on  the  minimization  of  (1)  considered  as  a  cost.  In 
order  to  avoid  getting  stuck  embedded  in  local  minima,  and  to 
accelerate  execution  time,  we  have  developed  a  feasible,  initial 
solution  to  the  optimization  scheme.  To  this  extent,  we  first 
compute  a  displacement  field  by  local  template  matching  on 
high  structure  points  and  subsequent  interpolation,  embbeding 
the  approach  into  a  gaussian  pyramid.  Our  implementation  re¬ 
lies  on  part  of  our  previous  work,  where  the  structure  detection 
is  based  on  the  eigenanalysis  of  the  correlation  matrices  of  the 
gradient  of  the  datasets  at  every  point  and  we  use  a  Kriging  [16] 
estimator  to  make  the  interpolation  [17],  [18].  It  is  interesting  to 
note  that  this  approach  is  able  to  deal  with  general  tensor  data, 
such  as  diffusion  tensor  MRI  and  hence  it  allows  the  generaliza¬ 
tion  of  the  method  proposed  here  to  tensor  data.  Nevertheless, 
for  the  purposes  of  this  paper,  any  template  matching  method 
which  is  able  to  provide  an  initial  solution  can  be  used. 

The  number  of  dimensions  of  a  3-d  displacement  field  search 
space  is  three  times  the  number  of  voxels.  In  order  to  save  com¬ 
putational  cost,  we  decouple  the  problem  searching  for  local 
minima  at  each  voxel.  The  optimization  algorithm  used  is  based 
on  a  Quasi-Newton  technique  called  BFGS  [19]  .  Separate  mini¬ 
mization  of  each  voxel  needs  knowledgment  of  neighbor  voxels, 
so  minimization  of  every  voxel  doesn’t  provide  us  the  nearest  lo¬ 
cal  minimum  in  the  global  search  space.  To  avoid  this,  several 
iterations  to  the  whole  dataset  are  carried  out  until  minimization 
reaches  a  consistent  solution.  To  avoid  false  optimizations  that 
the  regularization  term  cannot  deal,  a  median  filter  on  each  de¬ 
formation  field  component  is  carried  out  as  the  initial  field  for 
the  next  iteration. 

IV.  Two  Alternative  Methods 

In  what  follows  we  compare  the  proposed  method  with  the 
Euler-Lagrange  approach  and  with  a  Finite  Element  discretiza¬ 
tion.  A  brief  introduction  to  these  two  methods  is  given  below. 
Both  methods  are  implemented  using  a  gaussian  pyramid. 

1  Frobenuis  norm 


The  Euler-Lagrange  Method  consists  of  applying  variational 
calculus  to  equation  (1),  and  obtaining  a  list  of  differential  equa¬ 
tions  which  can  be  solved  iteratively.  This  procedure  is  a  diffi¬ 
cult  task  for  non-trivial  similarity  measures  and  regularization 
terms.  In  the  case  we  have  described,  the  similarity  measure 
and  the  regularization  term  have  been  approximated  to  first  and 
second  order  of  u. 

Using  Euler-Lagrange  equations  for  the  energy  functional  (1), 
we  reach  this  linear  system  of  equations,  defined  as 

ax  +  2(Anun+1  +  A12vn+1  +  A13wn+1)  =  2 a2(un  -  un+1)  (4) 
ay  +  2 (A21un+1  +  A22vn+1  +  A23wn+1)  =  2 a2(vn  -  vn+1)  (5) 
az  +  2(A31un+1  +  A32vn+1  +  A33wn+1 )  =  2 a2(wn  -  wn+1 )  (6) 

where  a\  and  A, ij  are  the  first  and  second  order  coefficients  of 
the  similarity  measure  respectively.  In  matrix  form  this  can  be 
depicted  as 

a  +  2Aun+1  =  2a2  (if1  -  iT+1)  (7) 

This  equation  can  be  solved  by  the  iterative  Gauss-Seidel 
Method  in  which  the  updated  components  of  the  deformation 
field  (wn+1,  vn+1 ,  wn+1)  are  computed  from  ( un ,  v n,  wn)  in  ev¬ 
ery  step  using 

w"+1  =  (A  -la2)-1^  -  la)  (8) 


The  Finite  Elements  method  [20]  relies  on  interpolating  the 
solution  using  as  nodal  points  the  vertices  of  a  discretization  of 
the  domain,  by  means  of  a  linear  formulation  where  the  values 
at  the  nodal  points  are  the  unknowns. 

In  the  3D  case,  we  have  selected  a  tetrahedral  mesh;  and  in 
the  2D  case,  we  have  chosen  a  triangular  mesh.  The  refinement 
of  the  mesh  is  chosen  as  a  function  of  the  level  of  detail  we  have 
assumed.  In  a  multi-resolution  scheme,  we  generate  a  mesh  with 
very  few  triangles  at  the  top  level  of  the  resolution  pyramid  and 
increase  the  number  of  triangles  as  we  move  to  higher  resolu¬ 
tions. 

In  order  to  accelerate  the  mesh  generation,  we  will  assign  big¬ 
ger  elements  to  less  detailed  zones  and  smaller  ones  to  high  in¬ 
formation  regions  of  the  dataset. 

The  algorithm  has  been  implemented  using  first  and  second 
order  approximations  of  the  correlation  coefficient.  Since  the 
results  derived  from  first  and  second  order  approximations  are 
not  significantly  different,  we  will  only  present  results  from  first 
order  approximations. 

V.  Results 

In  this  section  we  present  the  results  for  the  registration  of 
several  medical  datasets.  To  compare  the  performance,  we  have 
used  the  same  similarity  measure  and  the  same  regularization 
term. 

Figure  la  shows  a  small  region  of  an  axial  proton  density  (PD) 
MR  image  and  fig.  lb  shows  the  same  region  in  a  T2  weighted 
image,  though  slightly  deformed  by  a  synthetic  field.  The  syn¬ 
thetic  field  has  been  generated  by  smoothing  a  random  initial 
vector  field.  Figure  2a  shows  a  detail  of  the  synthetic  deforma¬ 
tion  field.  Figure  2b  shows  the  deformation  field  registered  by 
using  the  Euler-Lagrange  method,  with  a  =  0.8,  window  size 
w  =  5,  using  20  iterations  and  initialized  with  the  zero  solution. 
Figure  2c  shows  the  deformation  field  registered  by  using  the  fi¬ 
nite  elements  method,  with  a  =  0.8,  window  size  w  =  5,  and  a 
triangular  mesh  of  312  triangles  and  111  nodes.  Finally,  fig.  2d 


B.  Finite  Elements 


shows  the  deformation  field  registered  by  using  decoupled  op¬ 
timization,  where  the  initial  solution  has  been  obtained  using 
template  matching  with  a  search  window  of  size  wsearch  =  5 
and  a  similarity  window  of  size  wsimuarity  =  3.  The  local  op¬ 
timization  has  been  constrained  to  an  amplitude  of  five  pixels, 
a  =  0.8,  and  using  7  iterations  to  the  whole  dataset. 


a)  b) 

Fig.  1.  MRI  images  PD  (a)  and  T2  (b) 


c)  d) 

Fig.  2.  Synthetic  deformation  field  (a)  and  deformation  fields  obtained  with 
Finite  Elements  (b),  Euler-Lagrange  (c)  and  optimization(d) 


Figure  3a  shows  an  overlay  of  T1  w  and  T2w  MR  images  cor¬ 
responding  to  two  different  patients  and  the  surface  models  of 
their  ventricles.  Figure  3b  shows  the  estimated  deformation 
field  using  the  decoupled  optimization  approach  overlaid  onto 
the  models.  The  yellow  ventricle  is  considered  to  be  the  refer¬ 
ence,  and  the  green  ventricle  is  considered  to  be  the  current. 

In  order  to  compare  the  results  obtained  by  the  three  methods, 
fig.  4  shows  a  detail  of  one  coronal  section  of  the  reference  T1 
image  and  the  vector  fields  on  that  section  will  be  projected  onto 
it.  Figure  4b-d  shows  the  results  obtained  for  the  three  methods, 
EL  (a),  FE  (b)  and  DO  (c),  with  similar  parameters  to  those  used 
in  the  previous  experiment. 


Fig.  3.  Above:  ventricle  models  generated  with  the  MRI  T1  (green)  and  with 
the  MRI  T2  (orange)  from  two  different  patients  and  original  images.  Below: 
Three-dimensional  deformation  field  overlaid  on  a  ventricle  model.  Hot  colors 
means  large  deformations 


As  a  measure  of  quality,  the  absolute  correlation  coefficient 
has  been  computed  in  different  cases,  which  summarizes  table 
1 .  Results  for  a  synthetic  image  not  shown  are  also  displayed. 

VI.  Conclusions  and  Further  Research 

In  this  paper  we  have  introduced  a  novel  approach  to  general 
nonrigid  registration  problems  based  on  Decoupled  Optimiza¬ 
tion  (DO)  solution  to  a  variational  formulation  and  it  has  been 
compared  with  two  other  solutions  to  the  same  variational  prob¬ 
lem:  Euler-Lagrange  (EL)  and  Finite  Elements  (FE).  Although 
the  approaches  are  quite  different  in  nature,  they  all  give  simi¬ 
lar  and  meaningful  results  when  used  with  Correlation  Coeffi¬ 
cient  as  the  similarity  measure  and  eq.  (3)  as  the  regularization 
term.  The  main  differences  between  these  methods  are  related  to 


c)  d) 

Fig.  4.  T1  MRI  coronal  section  (a)  and  deformation  fields  obtained  with  Euler- 
Lagrange  (b),  Finite  Elements  (c)  and  direct  optimization  (d) 


TABLE  I 

Absolute  Correlation  Coefficients 


Registration 

Synthetic 

T2— PD 

Tl— T2 

Not  deformed 

1 

077S 

0.73  J 

EF 

0.84 

0.71 

0.69 

FE 

0.80 

0.70 

0.65 

DO 

0.85 

0.73 

0.70 

computational  cost,  and  possibility  of  incorporating  other  simi¬ 
larity  measures  or  regularization  terms.  The  computational  cost 
of  Euler-Lagrange  is  smaller  than  Finite  Elements  and  much 
smaller  than  Optimization.  Nevertheless,  the  DO  method  al¬ 
lows  an  easy  incorporation  of  any  similarity  measure,  such  as 
mutual  information  based  ones,  or  regularization  term.  More¬ 
over,  its  extension  to  deal  with  multivalued  data  is  also  possible, 
as  it  uses  algorithms  developed  to  deal  with  tensor  data  in  order 
to  get  good  initial  deformation  fields  and  the  optimization  solu¬ 
tion  is  powerful  enough  to  incorporate  the  tensor  reorientations. 
Nevertheless,  this  is  still  something  we  are  investigating  and  that 
will  be  reported  in  a  forthcoming  paper. 

Other  research  efforts  in  this  area  are  focused  on  comparing 
with  other  optimization  algorithms,  such  as  for  example  stochas¬ 
tic,  more  efficient  in  terms  of  computational  load  and  with  the 
ability  of  jumping  out  local  optima,  and  in  high  performance 
implementations  using  clusters  of  low  cost  personal  computers 
(beowulfs).  Our  motivation  is  to  obtain  a  robust  and  efficient 
computational  tool  which  can  be  easily  reconfigured  to  deal  with 
different  registration  problems  in  a  clinical  environment. 
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